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Abstract 

OEDIPUS is a Monte Carlo simulation program which can be used to determine the 
small- a; evolution of a heavy onium using Mueller's colour dipole formulation, giving the 
full distribution of dipolcs in rapidity and impact parameter. Routines are also provided 
which calculate onium-onium scattering amplitudes between individual pairs of onium 
configurations, making it possible to establish the contribution of multiple pomeron ex- 
change terms to onium-onium scattering (the unitarisation corrections). 
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PROGRAM SUMMARY 



Title of program: OEDIPUS 
Catalogue number: 

Program obtainable from: ftp: / /axpf . hep . phy . cam . ac . uk/pub/ theory/ oedipus . tar . gz , 
see also http : / / www . hep . phy . cam . ac . uk/ theory/software/ oedipus . html 

Licensing provisions: None 

Computers: Tested on Dec Alpha, Sun 

Operating system: Unix (OSF3.2, SunOS-4.1.3, Solaris-2.4) 

Program language used: Fortran-90 

Memory required to execute with typical data: 1-10 MB 

No. of bits in a word: 32/64 

No. of lines in distributed program, including test data, etc.: 6186 

Keywords: Small-x; BFKL; Pomeron; Onium-onium scattering; Dipoles; Monte Carlo; Uni- 
tarity. 

Nature of the physical problem 

BFKL evolution [1] of a heavy (quark-)onium, to give full information on the rapidity and 
transverse positions of gluons carrying a small fraction x of the onium's longitudinal momen- 
tum. This information can be used to calculate a variety of onium-onium scattering cross 
sections, including multiple pomeron contributions, which restore unitarity. 

Method of solution 

A Monte Carlo simulation of the evolution in rapidity of the gluon structure of the onium, 
using Mueller's colour dipole formulation of small-x physics [2]. 

Restrictions on the complexity of the problem 

The number of dipoles (gluons) in an onium grows as an exponential of the rapidity, and 
inversely with the value of a cutoff used to regulate an ultra-violet divergence. The time for 
evolution of an onium is proportional to the number of dipoles, while that for calculation of 
the onium-onium scattering amplitude goes as the square of the number of dipoles. Mem- 
ory restrictions also arise from the need to store large configurations of dipoles (very large 
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fluctuations can occur in the number of dipoles present). 
Typical running time 

On a DEC Alplia 3000 processor, evolution generates about 10000 gluons per second. At a 
rapidity y = 16 and an ultra-violet cutoff of 0.01 times the onium size, a day's running will 
generate about 10^ onium-onium scatterings. 
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LONG WRITE-UP 
1 Introduction 

There has been much interest recently in BFKL ||l|, |2|, |3| type processes and their possible 
observation at HERA and the Tevatron. Calculation of the cross sections for a number of 
these processes, such as diffractive dissociation (see e.g. [§, |5|) or exclusive vector meson 
production (e.g. p) requires detailed information about the dominant exchanged "object", 
known as the BFKL pomeron; for example one needs an understanding of the triple pomeron 
vertex and of the transverse distribution of small-x gluons. At present, the only BFKL type 
process which is fully calculable in perturbative QCD is high energy onium-onium scattering. 
Mueller's colour dipole formulation of this process |@, ^, |^, |l^ (for a related approach, see 
|11, |l2|), offers a well defined way of performing the necessary calculations, in the large Nq 



approximation, where Nq is the number of colours. 

The colour dipole formulation of small-x evolution and high energy onium-onium scat- 
tering is particularly suited to Monte Carlo simulation, since the evolution is probabilistic 
in nature, and because each branching (of one dipole into two) is independent of all other 
branchings. 

This paper is divided into three parts. The first gives a brief overview of the small- 
X dipole evolution of an onium and of the issues relevant to a Monte Carlo simulation of 
such an evolution. The second section examines the use of the dipole structures of a pair of 
evolved onia to calculate onium-onium scattering amplitudes, including the multiple pomeron 
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exchange contributions which restore unitarity. The final section gives a detailed description 
of the structure and use of the OEDIPUS package. 

2 Dipole structure of an onium 
2.1 Background 

One starts with an onium where the transverse separation between the quark and anti-quark 
is boi- The probability of generating a gluon at a point b2 carrying a small fraction of 
the light-cone momentum of the onium is : 



Nq is the number of colours {Nq = 3 for QCD), as is the strong coupling constant; A, the 
effective lifetime in rapidity (y) of the dipole, is related to the virtual corrections (or more 
intuitively, conservation of probability) : 



In the limit of large Nc, the colour structure means that the two new dipoles which arise 
(bo2 and bi2) can themselves independently emit gluons, while the original colour dipole is 
effectively destroyed. This branching of dipoles repeats itself until the rapidity of any new 
gluons which would be produced exceeds the maximum rapidity available. Through this 
branching process, one can establish the probability of any given configuration of dipoles. 

2.2 Notes on implementation 

One of the main aspects of eq. (||) is that it has non-integrable divergences at 602 = and 
612 = 0. This leads to the integral for being infinite. A solution to this problem is to 
introduce a cutoff on the dipole size, eliminating any region of b2 where 602 < p or 612 < p. 
This has to be used in both eq. (|l]) and eq. (|2|). The lifetime of a dipole of size h is therefore 
a function of h/ p. 

In the limit of small p the effects of the cutoff should disappear. However the number of 
small dipoles is large: the number of dipoles of size c from an onium of size 6, after evolution 
through y is approximately: 



dyd2b2 



dP 



(1) 




(2) 



(c,6,y)dlogc 



h 



exp((ap - l)y - \og{c/l)f /ky)d\ogc. 



(3) 



n 



CyJ-Kky 
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This is valid for |log(c/6)| <^ ky. The BFKL power is (a-p — 1) = 41og 2a5A''c/7r, and 
k = 14a5'A^c'C(3)/-7r, with C(3) ~ 1.202 being the Riemann zeta function. The smallest 
dipoles will be of size c ~ p, so as p is lowered the number of dipoles rises. In addition, the 
number of dipoles increases exponentially with rapidity. The time taken to generate an onium 
is proportional to the number of dipoles. The time to calculate onium-onium interactions is 
proportional to the product of the number of dipoles in the two onia. 

A final complication will arise because there are very large fluctuations in the num- 
bers of dipoles — the probability of obtaining a configuration with n dipoles goes as 
exp[— 7r(logn)^/4a5A'^C'y]) so that one has to allow for configurations which are very much 
larger than the mean. Since configurations generally need to be stored (for example to in- 
vestigate interactions between pairs of configurations), this can lead to considerable memory 
consumption, especially at large rapidities. 

The facility of imposing an upper limit on dipole sizes during the evolution is also in- 
cluded, to allow a crude investigation into the uncertainties due to infra-red effects. When 
the implementation of the upper cutoff is turned on, both dipoles produced from a branch- 
ing are required to be smaller than the upper cutoff. The lifetimes of dipoles are adjusted 
accordingly. 

The evolution is carried out with fixed as, which is justified in the limit of very heavy onia. 
There is no known unique way of including a running coupling constant, however it would be 
possible to modify the program to implement a scheme such as that used by Nikolaev and 
Zakharov 



3 Onium-onium scattering 
3.1 Background 

One wishes to obtain the amplitude F{h,h' ,r,Y), for elastic scattering between two onia of 
sizes (and orientations) b and b', whose centres are separated by a transverse distance r, with 
a total rapidity between them of Y (~ logs, where -y/s is the centre of mass energy). 

Let 7 be a particular dipole configuration for an onium, which contains dipoles of position 
and size (ri, ci) . . . (r„^, c„^). The interaction between two such onia, moving in opposite 
directions is ^, |l^ : 

A.y = im/(i^i-i^j'Ci,c^-)> (4) 

1=1 j=l 

where /(r, Cj, c' ), the interaction between dipoles of size c and c', whose centres are separated 
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by r is I, 



/(r,c,c') = ^ 



al r |r + c/2-c72||r-c/2 + cV2| 
log 



n 2 

(5) 



|r + c/2 + c72||r - c/2 - cV2 
One has to average f^^y over dipole configurations 7, 7' of the two onia. Therefore at the 
level of one pomeron exchange, the elastic amplitude is: 

F(i) (b, b', r, y) = - ^ P^(ro, b, y)Py (ro + r, b', Y - y)f^^y, (6) 

7:7' 

where P^(ro,b,y) is the probability of obtaining a dipole configuration 7 after evolution 
through y of an onium of size b centred at tq. Note that in this case the result can be shown 
to be independent of y, the division of rapidity between the two onia (or equivalently of the 
longitudinal frame in which the calculation is performed). 

One can also calculate contributions to the amplitude involving the exchange of k pomerons 



F« (b, b', r,Y) = ^ Y. ^7(ro, b, y)Py (ro + r, b', Y - y)i-f^,y)K (7) 

These multi-pomeron terms, though formally non-leading in 1/Nc, are enhanced at high 
rapidities by the large numbers of dipoles in each onium. All numbers of pomeron exchange 
can be resummed to give an amplitude F which explicitly satisfies the unitarity bound: 



(b, b', r, y) = -1 ^ F,(ro, b, y)Py (ro + r, b', Y - y){l - e-^.V ). (8) 

These last two equations are approximations based on there being a large number of dipoles 
interacting. Multiple pomeron contributions should be calculated in the frame y = Y/2, where 
corrections that would arise from wave-function saturation, which is not being calculated, are 



expected to be smallest (see also [15]). 



3.2 Implementation considerations 

One of the main difficulties in calculating onium-onium scattering amplitudes is the large 
numbers of interactions which have to be worked out: in principle the product of the numbers 
of dipoles in each onium, multiplied by the number of points at which one needs to know the 
interaction. 

The first way to ease the problem is to note that the dipole-dipole interaction, eq. (|5|), is 
relatively local: at large distances it dies off as where r is the separation between the 
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centres of the dipoles. Therefore it should be safe to neglect the interaction between dipoles 
if their separation is sufficiently large, say greater than twice the sum of their sizes. This is 
found to reproduce the total dipole-dipole interaction to better than one percent^. 

One generally wants to know the amplitude at many different impact parameters (i.e. 
onium-onium separations), for example to work out the total cross section (with the normal- 
isation used above, just twice the amplitude, integrated over all impact parameters). One 
way to do this is to work out the interaction on a grid. Since there are large fluctuations in 
the spatial extent of the dipole configurations, the grid size has to be variable — it should 
completely cover the area where dipole-dipole interactions are significant. By summing the 
interactions at each point of the grid one obtains an approximation to the total amplitude 
for that configuration pair (the finer the grid, the better the approximation). 

This approach is necessary if one wants to find a good approximation to the total amplitude 
for each configuration pair, however it is slow and does not offer any easy way to store the 
results so that they can be retrieved for later analysis (say if one only wanted to look at 1 
and 2 pomeron exchange, but after the run, decided that 3 pomeron exchange would also be 
interesting) — if one uses a grid which is 100 x 100, then one requires 40 kilobytes per pair 
of configurations, for a single value of rapidity. To gain adequate accuracy one needs several 
tens of thousands of events; with, say, 5 values of rapidity this gives several gigabytes of data. 

One reason why it is important to sample a large number of configuration pairs, is that 
rare, large configurations contribute significantly to the total cross section 0. Multiple 
pomeron contributions are dominated by rare, dense configurations. For a given pair of 
dipole configurations, however, the variation of the amplitude (at those points where it is 
significant) with impact parameter is not so large. So it can be advantageous to evaluate 
fewer points in impact parameter (i.e. making a worse approximation to the amplitudes of 
individual configuration pairs) in exchange for a larger sample of configurations. The limit 
of this is to evaluate the onium-onium interaction, not on a grid, but along a radial line. In 
addition, because there are fewer points in impact parameter it becomes tractable to store the 
results. In fact the individual amplitudes are not stored event by event, rather a probability 
distribution for the amplitude is stored at consecutive bins in radial position. This information 
can then be retrieved and used to determine whichever quantity one is interested in. 



^Note that this may not correctly reproduce the profile of the interaction at very large onium-onium 
separations, where the density of dipoles is also dying off also as 1 /r'* (though enhanced by double leading 
logarithmic factors). However this region is not important for the total interaction. 
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4 Structure and use of OEDIPUS 



The program comes in many different files in various directories: basic_src/ contains all 
the common subroutines which are involved with data input and output, initialisation and 
evolution of onia, as well as routines for determining onium-onium interaction amplitudes. It 
does not contain any code for analysing the results from the evolution or interaction: this 
must be provided as a "driver" program by the user. A number of examples are to be found 
in the directory samples/, including drivers to examine the spatial distribution of dipoles in 
the onium wave function and to determine the total and elastic scattering cross sections. 

4.1 Storage of gluons and dipoles and evolution 

Gluons and dipoles each have a specific data type: 

type gluon 

complex (kind ( IdO) ) :: psn ! stores the transverse position 

double precision : : rapidity 

end type gluon 

Complex variables are used throughout the program to deal with transverse positions. The 
dipole type is 

type dipole 

type (gluon), pointer :: lo_y_gluon, hi_y_gluon 
double precision : : size 

type (dipole) , pointer :: child_l, child_2 
integer : : index 

end type dipole 

It is useful to separate the gluon content of the dipole into high rapidity and low rapidity, 
mainly for ease of coding. In the case of the initial onium, the two quarks are represented by 
gluon types. The child_l and child_2 pointers are needed to store the tree structure of the 
evolved onium. If they are not associated then the dipole has no descendents (for example, 
because any descendents would have a rapidity larger than the maximum of the evolution). 
Each dipole in an onium has a different index — this can be useful in analysing the branching 
structure. 

The following sequence of routines is needed to produce an onium and evolve it, for each 
loop (or "event") of the program. First, various internal counters must be reset by calling 
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subroutine reset_counters 

Then an onium is produced with 

subroutine init_onium(oniuni, onin_size, phi) 

This returns a pointer onium to a dipole of size omii_size, centred at (0,0), with orientation 
phi. To produce a randomly oriented dipole, call 

subroutine init_onium_rnd_phi (onium, onm_size) 

The onium is evolved up to a rapidity maxy by calling 

subroutine evolve_onium(onium, meixy, n_tree) 

where, on return, n_tree is the total number of dipoles in the tree. To be able to analyse 
the results, one wants only the dipoles associated with a particular rapidity y, so one must 
extract them with: 

subroutine extract_onium_dpls (onium, y, extrctd_dpls , n_extrctd) 

The array extrctd_dpls ( : ) is of type dpl_pntr (arrays of pointers are not permitted in 
Fortran 90 — the solution to this is to define a type which contains nothing but a pointer), 
and on return contains pointers to the n_extrctd dipoles. The array must be sufficiently 
large: one can ensure that it is, by dynamically allocating it to be of size (n_tree / 2 + 1). 
The array of dipole pointers can then be processed by the Tiscr. 

This sequence allows one to first evolve the onium up to the maximum rapidity that one 
is interested in, and then examine the dipole structure over a range of rapidities. 

Additional onia (for use simultaneously) can be produced, evolved, and the dipoles ex- 
tracted, by repeating the same set of calls (omitting the call to reset_counters). 

4.2 Onium-onium interaction 

Two sets of routines arc provided which calculate amplitudes for onium-onium scattering. 
The first returns a reasonable estimate for the elastic scattering amplitude between a pair of 
onium configurations, for all relevant impact parameters: 

subroutine onm_onm_grid(ext_dplsl , nl, ext_dpls2, n2, grid, gh_lo, & 
& gh_int, gv_lo, gv_int) 
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It takes a pair of dipole sets (nl dipoles in ext_dplsl and n2 dipoles in ext_dpls2) and 
determines the range of relative separations over which they interact. The onium-onium 
interaction is then determined for each point of a grid covering this area (grid(: , :) must 
be a two-dimensional square array; its size determines the resolution of the sampling) — the 
program returns information relating the array to the physical grid being used: the horizontal 
and vertical positions of the bottom left hand edge of the grid (gh_lo and gv_lo), as well as 
the horizontal and vertical spacing (or interval) of the grid (gh_int and gv_int). For each 
dipole pair, it only works out the contribution to those points on the grid where the dipole- 
dipole separation is not too large (as discussed in section 3.2). The amplitudes are always 
worked out at the centres of the grid rectangles. 

The second approach works out a binned probability distribution for the amplitude in 
many sections along a radial line. The details of the binning are held in the binning module 
which provides a variable bn of type bin.prms: 

type bin_prins 

! binning vars for amplitude prob distributions 
integer :: n_f_lo, n_f_hi, nf 

double precision : : f_inid, f_lo_lgint , f_hi_lgint 
! binning vars for r 

integer : : n_r_lo, n_r_hi, nr 

double precision : : r_mid, r_hi_lgint 
end type bin_priiis 

The binning of the amplitude probability distribution is divided into two regions: the lower re- 
gion is for amplitudes below f jnid and the bins are indexed to (n_f _lo - 1), logarithmically 
spaced with an interval f _lo_lgint. It is important to have a bin explicitly for the absence 
of interaction (otherwise when integrating the interaction over a large area the smallest bin 
amplitude could contribute significantly) — this is bin —1. The lowest non-zero bin must be 
sufficiently small to accurately reproduce the integral of the dipole-dipole interaction (bearing 
in mind that because the interaction dies off as there are regions where the amplitude 

can be quite small, but whose overall contribution to the amplitude is non-negligible). It is 
to allow adequate coverage of this wide range of small amplitudes that logarithmic binning 
(in the small amplitude region) is used. For amplitudes above f _mid the binning is performed 
logarithmically with an interval f Jii_lgint. The large amplitude region will contribute more 
significantly to the total amplitude and dividing the binning into two regions allows one to 
sample it more accurately in the regions which are more important. The reason to have loga- 
rithmic binning for large amplitudes is to be able to store the very large amplitudes which can 
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occur at large rapidities. It is sometimes useful to have the index of the highest amplitude 
bin, which is nf = n_f _lo + n_f Jii. The variable bn is initialised with default values which 
reproduce total amplitudes to an accuracy of about one percent. 

The binning in r, the distance along the radius, is also divided into two regions, but for 
different reasons: there are occasional events with very large transverse extents. To be able 
to include them with linearly spaced bins would require a prohibitively large number of bins. 
So for large r, logarithmic spacing of the bins is useful. But for small r, the most efficient 
sampling is one where each bin corresponds to a region of similar area — linear binning in r 
is more appropriate there. There are bins numbered to (n_r_lo - 1) linearly spaced from 
r = to rjnid. Above this, there are n_r Jii logarithmically spaced bins, with an interval 
r_lii_lgint. 

To help bin and "unbin" the data, the following functions are provided: 
function bin_of_f (f) 

returns the bin number associated with an amplitude f . If the amplitude falls into a bin 
higher than nf then it is put into bin nf — this serves as a form of overflow procedure: if the 
highest bin contains any entries then the program should have been run with a larger number 
of amplitude bins. The inverse function 

function f _of _bin(bin) 

returns the value of the amplitude corresponding to the (logarithmic) centre of the bin (on 
average, the optimal value to use in reconstructing amplitudes). 

For a given dipole-dipole pair, it is useful to know to which radial bins their interaction 
will contribute. The function 

function irad(r) 

aids this by returning the bin corresponding to a given radius (the interaction region for a 
pair of dipoles is defined by two radii, which are converted into the range of bins which must 
be sampled). For each radial bin, it is useful to sample at a random radius within the bin. 
The random value of r returned by 

function rad_rnd(i) 

has a probability distribution which increases linearly with r within the bin i, replicating the 
distribution of radii that would be obtained by choosing a random position in the correspond- 
ing ring. Finally it can be necessary to know the lower radius of a given bin (for example 
when analysing binned results), which is obtained by calling 
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rad_lo(i) 



Given a pair of dipoles sets from evolved onia (nl dipoles in ext_dplsl and n2 dipoles in 
ext_dpls2), the following routine 

subroutine add_f _to_bins (ext_dplsl , nl, ext_dpls2, n2, fbins_single) 
adds 1 to the appropriate amplitude bin of 
f bins_single (-1 :bn°/onf , :bn7onr) 

for a radius chosen randomly with rad_rnd in each radial bin. 
4.3 Data input and output, and initialisation 

A wide variety of information can be extracted from the evolved onium and the program 
needs to have the flexibility to input and output whatever information the user wants. In 
addition, it must be possible to restart a Monte Carlo run to add to previously determined 
data — this should be as automatic as possible. A number of routines are provided to make 
this simpler. 

read_parains_f dat (onm_size , maxy, n_prev_events , n_new_events) 

This routine first parses the command line. Defining argn to be the n*'' command line argu- 
ment, it opens argl .dat (a file containing formatted data from the evolution) and argl .prm 
(a file containing the parameters of the evolution). The number of evolutions (or events) to 
be performed in this run of the program (n_new_events) is the value of arg2. It then looks 
to see if arg3 is present: if not (or if it is a dash) then the evolution is a continuation of a 
previous one and the contents of argl .prm are read in, in the following format: 

jseed(l) jseed(2) 
r_cut_lo r_cut_hi 
onm_size maxy 
data_io 
n_prev_events 

with the following definitions: 

jseed(l:2) the pair of integers used as a seed for the random number generator HWRGEN, 



which uses I'Ecuyer's method, described in as implemented in the HERWIG ||17| 
event generator; 
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r_cut_lo the lower cutoff on dipole sizes; 

r_cut Jii tfie upper cutoff (if it is negative, tfien no upper cutoff is implemented); 
onm.size the onium size; 

maxy the maximum rapidity up to which onia are evolved; 

data_io a variable provided by the user in the module data_prms holding any additional 
information needed by the user's program for treating the data (e.g. a range of dipole 
sizes being studied). It can be of any type (including user defined type); 

n_prev_events the number of "events" (i.e. Monte Carlo evolutions) processed so far. 

The routine also performs the initialisation which is necessary before evolution can be carried 
out: a call to init_lives which sets up an array tabulating dipole lifetimes as a function 
of their size, and a call to nullif y_child_set which nullifies pointers to child dipoles (since 
Fortran 90 has no mechanism for specifying the initial status of a pointer). It is the user's 
responsibility to read in the data from unit idat (since the format and nature of the data 
will vary from program to program). 

If argS is present then the file starters/arg3 (path from top directory of the distribution) 
is read in. The format is the same as the first three lines of argl.prm. If either value of 
jseed(l:2) is zero then the default seed is used; n_prev_events is set to zero. It is then the 
user's responsibility to initialise any data variables and provide default values for data_io. 
This mechanism makes it much easier to run the same program with different values of cutoffs 
and rapidity (the most commonly varied parameter) by having a single set of starter files in 
the the starters/ directory — effectively by a command line argument, rather than by 
having to edit a new parameter file each time. A number of starter files are provided with 
the distribution and are described in the file starters/README. 

write_params_f dat (onm_size , maxy, n_events, param_arg, dat_arg) 

This routine is called to output the updated parameters of the evolution (usually only the 
seed and the number of evolutions so far, n_events, will have changed). It rewinds both 
the parameter and data devices. Normally the parameter and data devices (paraiii_arg and 
dat_arg) will be the same as those used for input (iprm and idat) but the arguments are 
available (though optional) because there is occasionally a need to use different output devices 
(e.g. if the user wants to retain the original data). It is the user's responsibility to output the 
data. 
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The advantage of using formatted data is that it is easy to examine or plot. But it is 
very inefficient for storing large amounts of data. Two routines are provided for dealing with 
unformatted data (the parameters file is still formatted): 

read_params_uf dt (oiim_size , maxy, n_prev_events , n_new_events) 
write_params_uf dt (onm_size , maxy, n_events, param_dev, ufdt_dev) 

They are very similar to the previous routines. The main difference is that the data device 
is now iufdt, which is opened for unformatted output. Since the most common use of 
unformatted input and output is for the binned probability distribution of the amplitude 
(see section 4.2) the parameter file contains the additional set of parameters: bn (a variable 



in the binning module, of type bin_prms) which defines the spacing of the bins. With the 
unformatted data files it will often be necessary to read them in again to analyse them. 
Because the various parameters may be needed to know how to read in and process the data 
a routine 

read_params_old_uf dt (prnt_size , maxy, n_events) 

is provided. It looks at argl to determine the files to use and opens them (with the OLD 
specifier, to prevent them being accidentally altered). It does not perform any initialisation 
for the evolution. 

For binned amplitude data, the format is fairly constant, so routines are also provided for 
reading and writing the data: 

subroutine allc_read_f bins (n_y , n_prev_events) 

This allocates an array 

fbins(-l :bny„nf , 0:bnyonr, n_y) 



in accord with the format for binned data discussed in section |4.2| reads it in from device 
iufdt. The extra dimension comes in because one might want to store results for several 
(n_y) values of rapidity from a single run. The reason for the argument n_prev_events is 
to produce the correct normalisation for the data: on disk, it is stored so that each bin 
contains the probability associated with it. When adding extra data, it is convenient to 
use a normalisation where each bin contains the number of events contributing to it. If 
n_prev_events is then the bins are not read in, but instead initialised to zero. If the 
argument is not present then the normalisation of the bins is not changed. 
To write the data, use the routine 
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subroutine write_f bins (n_events , iufdt_arg) 

It is sometimes useful to be able to specify a different device, iuf dt_arg (optional), for output 
than for input, for example when writing out only a portion of the data. The number of events 
needs to be provided so that the stored normalisation on output is consistent with each bin 
containing a probability. 

The data is stored as real rather than double precision to reduce the use of disk space. 



4.4 Some internal details 

The gluons and dipoles which are generated by evolution have to be stored, for example so 
that one can then work out onium-onium interactions. One of the main problems is that one 
doesn't know beforehand how many dipoles will be produced — the mean number of dipoles 



can be estimated, however the fluctuations above this are very large (see section 2.2). A 
solution is to dynamically allocate memory for gluons and dipoles as the evolution proceeds. 
This slows down the evolution (by nearly a factor of two on some systems) because of the 
time needed for allocation and deallocation. 

The alternative is to provide fixed size arrays for the dipoles and gluons, where the size is 
chosen beforehand — essentially by trial and error. For cases where the largest configurations 
will be using a significant fraction of a machine's resources, this has the disadvantage that 
by having previously allocated the memory, it is taken up for the whole duration of the job 
(even if there is nothing stored there), whereas with dynamic allocation, the memory is used 
up only when it is needed. 

Limits have been coded into the gluon module for the maximum number of gluons: the 
variable inax_gluons. For the situation with dynamic allocation, this is the maximum number 
of gluons in each onium. It is present simply to allow the program to handle cases which might 
otherwise cause it to crash because they used all the available data store for that process. In 
the case of the fixed arrays, the quantity represents the maximum number of gluons in total 
(i.e. from all evolved onia). Apart from that, when using the subroutines described above for 
initialising onia, evolving them, etc., the functionality should not depend at all on the method 
of gluon storage. 

There are ways round the limitations on the maximum number of dipoles. One could store 
very large configurations on disk (the computer would normally be doing this anyway with 
the swap space — but if the program is responsible then the swap space is still left for other 
processes). Alternatively, there are situations where the dipole structure need not be stored 
at all. The occasions when a stored dipole structure is needed is when looking at some kind of 
correlation between dipole pairs (such as onium-onium interaction), because for each dipole, 
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one needs access to all the other dipoles. But in these situations, the time per configuration is 
normally the limiting factor, since this is proportional to the square of the number of dipoles. 
In cases where the time taken for some analysis is linear in the number of dipoles, often, it 
is not necessary to have knowledge of the whole configuration. So an evolution routine could 
be written which called an analysis routine each time it reached the "tip" of a "branch" of a 
dipole "tree" , and it would never have to store more than a single path from the base of the 
tree to the furthest "tip" . 

4.5 Compilation 

The program is provided in the form of a number of separate files, the contents of which 
are documented in README files. To aid in the compilation process, two shell scripts are 
provided which select the correct set of source files. Detailed information on their functioning 
is provided as part of the distribution, however a summary is given here: 

mkf driver [ dyn | sml | Irg ] 

This is for compiling a driver routine which uses formatted data i/o. The name of the 
driver should be provided without the .f90 extension. The second command line argument 
indicates the kind of gluon and dipole storage used — the default is to use dynamic allocation 
(corresponding to the option dyn). Two fixed size array storage options are also provided: sml 
(up to 5000 gluons) and Irg (up to 400000 gluons, which with a lower cutoff of 0.016, should 
be adequate for evolution up to rapidities of y ~ 15). The executable is named driver_dyn 
(or sml or Irg as appropriate) to allow multiple copies of the executable, with different gluon 
storage, to coexist. 

mku driver [ dyn | sml | Irg | noev ] 

This is for compilation of drivers which use unformatted data i/o. It also compiles in the 
routines and modules associated with binning of the amplitude. The command line arguments 
are the same as for mkf, except for the extra option available as the second argument: noev. 
This is to be used for drivers which analyse data from a previous evolution — no evolution 
routines are compiled, and the data and parameter files are opened with the ' OLD ' specifier, 
meaning that they cannot be modified. Naming of the executable follows the same convention 
as for mkf except when the noev option is used, in which case the executable is named driver 
(since there should be no need to have multiple copies compiled with different gluon storage 
methods) . 
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4.6 Test run 



The run described here tests the evolution of the onia, determination and binning of the 
interaction amphtude, and subsequent analysis of data for onium-onium interactions. It is for 
moderate rapidities and very small statistics to ensure that it runs quickly (a few seconds) 
on most systems. It is necessary to first compile the evolution routine (to be found in the 
sample/ directory of the distribution) with 

> mku onm-onm-bnd dyn 

which will include in dynamic allocation of memory for gluons and dipoles. The first command 
starts an evolution (with 10 events) where the maximum rapidity for each onium is 4, giving 
a total rapidity for the collisions of 8. The lower cutoff used is 0.1 times the onium size. The 
binned, unformatted data are stored in test_y8.ufdt (beware: this file is over 1MB in size) 
and the parameters in test_y8 .prm. The format of the variable data_io (together with the 
defaults used here) is 

data_ioy„y_wvf n = OdO 
data_ioyoy_int = 0.5d0 
data_ioyon_y = 5 

which specifies that the onium-onium interaction is determined for 5 values of rapidity, each 
separated by 0.5 (note that this is for each onium — so the intervals in total rapidity are 
1.0). The maximum rapidity at which dipoles are extracted from each onium is maxy, but 
each onium is evolved to data_io7„y_wvf n (or maxy, if this is larger) allowing different runs 
to effectively use the same dipole configurations but to examine the interactions at different 
values of rapidity. Various messages are output at the start of the program, as different parts 
of the startup procedure are accomplished. The line dealing with backups may be preceded 
by other messages on some systems (e.g. those using the NAG compiler — see section [4.8[ ), 
because the compiler does not offer the facility of running a system command. 

The next command causes the evolution to be continued for a further 5 "events" (simply 
to allow a test of the continuation function). One then wants to analyse the results of the 
evolution. In this test we will look at the total amplitudes (the onium-onium amplitude 
integrated over all impact parameters, which is equal to half the onium-onium total cross- 
section). The program to do this is compiled with 

mku ftot noev 
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The option noev ensures that unnecessary routines (such as those required for evolution) 
are not compiled in. The output includes various messages at the beginning indicating the 
initialisation procedure (these messages are sent to standard error, while the numerical output 
is sent to standard output, to allow it to be redirected to a file). The format of the numerical 
output from noev is one line for each value of rapidity, where each line has: 

(total rapidity) (unitarised amplitude) (1 pomeron amplitude) . . . 
(n_pom pomeron amplitude) 

(all on one line), and n_pom = 4 is a parameter in the program. It is the absolute values of 
the amplitudes which are output. 

A second test run, which tests the formatted data output is provided and documented 
within the distribution. In addition, some other driver and analysis routines are provided. 
These can all be found in the samples/ directory. 



Examples of results produced with OEDIPUS can be found in [13, 14]. 
Typical rapidities for which adequate statistics are accessible with a day's running are 
Y < 20, using a lower cutoff of 0.016. 



4.7 Distribution 

The distribution (provided as a compressed tarred file) consists of a number of directories. 
Details of their contents are documented in extensive README files. A shell script is provided 
to help automate the installation procedure. Instructions are also included for installation by 
hand. 



4.8 Portability 

The main issues of portability relate to accessing command line arguments and causing a shell 
script (which performs backups of data and parameter files) to be executed. Interfaces to the 
system routines to perform these functions are provided as routines whose names start with 
the lcl_ prefix, held in a file 

basic_src/ common/lcl_??? . f 90 

where ??? may be dec, nag or std (or the user can generate a new file for his/her operating 
system). The file to be used must be specified as part of the installation procedure. With the 
NAG compiler (version 2.1), there does not seem to be any way of causing a shell script to be 
executed (the usual unix system subroutine seems to be absent). The routines for accessing 
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command line arguments while available, do not seem to be documented. If the Fortran 90 
compiler offers no mechanism for accessing the command line arguments then one can use 
lcl_std.f90, which prompts the user for the information held in each of the arguments. It 
conforms to the standard. 

Local implementations are also necessary for writing large, unformatted arrays to disk (on 
some systems, the i/o buffer is not sufficiently large and the array must be written in small 
sections). Details are provided in the lcl_???.f90 files. 

To port to a non unix system (such as VMS), it would in addition be necessary to rewrite 
the compilation scripts. 

4.9 Global data 

All global data is stored in the form of modules. These are summarised in table |l] (numbers 
in brackets refer to default values). 

In addition there are modules associated with the binning routines. These are summarised in 
table ll 
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Module name 


V C\i± lC\il_f 


T^osrn T~if 1 on 


constants 


a ~\ TiVl a C! 


f-X-y 1 tliC- OLiVJli^ (^IJ U.]Jilil^ L^Uilo UdliU lU.Xi III ' ' ') 
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fa ptor of AT^-t f — a 1 nh a q i 


cuts 


T* <"n^ ~\ n 
■L U. U Lvy 
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UllC lU W Cl U. t 'Jll 


T* r'li'l" ^^ 1 

■L U. U _JJ. J- 


tnP IITITIPT* PlltOTT 
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hi_cut_impleiii 


true if an upper cutoff is being used 


J. (J 111 U LL U. J- 


J. UJ. ill 


1 / o riPAri/^P Tdr T»tiT*£imptPT'Q i 1 D 1 
1 / u VJ.C V ii^c lui ^di ojiiic tci o I X u J 


id.cL't 


i/o Hpvipp for formattpH ontont M 2^ 


iufdt 


i/o device for unformatted data (14) 


gluons 


type gluon 


the type definition for gluons 


type dipole 


the type definition for dipoles 


type dpl_pntr 


type definition for a pointer to a dipole 


max_gluons 


the maximum number of gluons to be generated 


data_def s 


data_io 


variable of user chosen type, with details of data i/o. 
This module must be provided by the user. 


sub_def s 




interface blocks for all the main subroutines 



Table 1: A list of the modules (and their main contents) needed by user's programs 



binning 


type bin_prms 


type definition to store all the parameters for binning 
of amplitudes 


bn 


variable actually containing the information 


f binsjndl 


fbins 


array containing amplitude bins 


bin_interf ace 




interface blocks for the binning routines 



Table 2: Modules associated with routines for binning the amplitude probability distribution 
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TEST RUN 

The input and output of the test run are shown below: 



> onm-onm-bnd_dyn test_y8 10 y4 
Starting a new evolution 
HWRGEN will use its default seed 
Lifetimes have been initialised 
Have allocated bins 

If there are no error messages then the backups have been done 
Am about to evolve 

> onm-onm-bnd_dyn test_y8 5 
Continuing a previous evolution. . . . 
HWRGEN has been initialised 
Lifetimes have been initialised 
Have allocated bins 

If there are no error messages then the backups have been done 
Am about to evolve 

> ftot test_y8 

Reading from a previous evolution. . . . 
HWRGEN has been initialised 
Have allocated bins 

4.0000E+00 1.3183E-01 1.3466E-01 2.9205E-03 9.6577E-05 3.3521E-06 
5.0000E+00 1.8955E-01 2.0176E-01 1.4294E-02 2.4481E-03 4.1440E-04 
6.0000E+00 2.0235E-01 2.1490E-01 1.4465E-02 2.2114E-03 3.4580E-04 
7.0000E+00 4.2093E-01 4.7515E-01 6.3966E-02 1.1441E-02 1.9750E-03 
8.0000E+00 4.4738E-01 5.0385E-01 6.5855E-02 1.0901E-02 1.7452E-03 

> 

The starter routine y4 is the following (default seed, lower cutoff of 0.1, no upper cutoff 
(—2.0), onium size of 1.0, and maximum rapidity for each onium of 4.0): 


0.1 -2.0 
1.0 4.0 
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